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SUMMARY 

An  adaptation  of  the  two-step  Lax-Wendroff  method  is  used  for  solving  the 
unsteady  one-dimensional  equations  of  non-linear  j hallow  water  theory,  including 
both  frictional  resistance  and  lateral  inflow  terms.  This  finite  difference 
method  is  fast,  accurate  and  simple  to  programme  and  covers  the  formation  and 
subsequent  histoiy  of  discontinuities  in  the  solution,  in  the  form  of  bores  and 
hydraulic  jumps,  without  any  special  procedures.  The  behaviour  of  the  numerical 
solution  behind  these  jumps  is  found  in  the  examples  to  be  sufficiently  smooth 
without  the  addition  of  an  artificial  viscous  force  term.  A  variety  of 
illustrative  examples  is  given,  including  simple  cases  of  flood  waves  in  rivers, 
bores  in  channels  resulting  from  rapid  changes  of  upstream  conditions, 
oscillatory  waves  cn  a  super-critical  stream  and  a  simple  hydrology  example  with 
a  significant  lateral  inflow  from  rain.  Several  checks  of  the  numerical  method 
are  included.  The  examples  are  confined  to  channels  of  uniform  rectangular 
cross-section,  but  the  method  generalises  in  a  straightforward  way  to  real  rivers 
end  estuaries  in  which  the  cross-section  is  non-re ct angular  and  varies  along  the 
length  of  the  channel. 
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Many  practical  problems  of  open  channel  hydraulics  can  be  modelled  by  the 
equations  of  unsteady  one-dimensional  non-linear  shallow  water  theory, 
particularly  if  a  frictional  resistance  term  (e.g.  Ch^ssy)  and  a  lateral  inflow 
term  are  included.  For  example,  flood  waves  in  rivers,  surges  travelling 
along  channels,  tidal  flooding  in  estuaries,  and  nurfact  run-off  from  heavy 
rain.  Shallow  water  theory  is  applicable  to  flows  in  which  the  wave  length  of 
disturbances  and  the  radius  ^.f  curvature  of  the  water  surface  is  much  greater 
than  the  depth  of  water;  •non-linear*  means  no  restriction  is  placed  on  the 
ratio  of  wavs  amplitude  to  water  depth. 

The  equations  of  shallow  water  theory  are  similar  to  those  of  gas  dynamics. 
The  purpose  of  this  work  is  to  adapt  a  method  that  has  been  found  very  success¬ 
ful  in  the  gas  dynamics  application,  namely  the  two-step  Lax-Wendroff  method 
(Richtmyer  and  Morton1),  to  the  shallow  water  equations.  The  main  differences 
are  that  the  shallow  water  equations  are  not  in  1 conservation-law 1  form  (e.g. 
df/dt  «  dg/dx)  whan  frictional  resistance  and  lateral  inflow  terms  are  included 

m  tm 

and  that  there  are  just  two  equations  (continuity  and  momentum)  against  the 
three  of  gas  dynamics  (continuity,  momentum  and  energy) .  In  shallow  water  flows, 
bores  and  hydraulic  jumps  (stationary  bores)  correspond  to  shock  waves  in  gas 
dynamics.  These  points  of  discontinuity  of  the  mathematical  solution  are 
referred  to  collectively  as  •jumps*. 

The  equations  of  these  flows  are  far  too  complicated  for  an  analytical 
solution  to  be  possible  and  a  numerical  method  is  essential.  The  first  choice  is 
between  an  Eulerian  and  a  Lagrangian  formulation.  The  latter  has  advantages  in 
gas  dynamics  when  mixtvres  of  gases  with  differing  thermodynamic  properties  are 
involved,  but  the  Eulerian  form  ia  more  convenient  in  the  current  application; 
•also  an  Eulerian  method  is  more  easily  extinded  to  unsteady  two-dimensional 
problems.  The  second  choice  is  between  a  characteristic  finite  difference  method 
and  a  direct  finite  difference  method.  In  the  latter  methods  the  continuity 
and  momentum  equations  are  expressed  directly  in  terms  of  finite  difference 
approximations,  while  in  the  former  it  is  the  characteristic  properties  of  these 
equations  that  are  so  approximated.  Numerical  solutions  based  on  characteristics 
are  accurate  and  unconditionally  stable  numerically,  and  thinking  in  terms  of 
characteristics  helps  both  physical  interpretation  and  seeing  whether  a  problem 
with  given  initial  and  boundary  conditions  is  ‘well-posed*.  But  such  methods 
are  slow  and  become  complicated  when  the  flow  contains  jumps,  particularly  when 
a  jump  is  moving  into  fluid  that  is  already  disturbed.  On  the  other  hand,  finite 
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difference  methods  are  fast  but-  until  the  Lax-Wendroff  method,  were 
relatively  inaccurate.  The  Lax-Wendroff  method  has  second  order  accuracy  and 
has  been  found  in  gas  dynamics  to  be  applicable  to  flowe  containing  Jumps, 
without  special  procedures,  such  as  shock  fitting  or  ucing  an  artificial  viscous 
force,  being  essential  as  with  other  methods  to  avoid  violent  instability  in 
the  numerical  solution  behind  a  Jump.  The  Jump3  appear  automatically  when 
using  the  Lax-Wendroff  method  on  these  flows,  as  near-discontinuities  across 
which  the  dependent  variables  have  very  nearly  the  correct  Jump  and  which 
travel  at  very  nearly  the  correct  speed  through  the  fluid.  The  jump  is  spread 
over  about  four  finite  difference  space  steps,  with  a  slight  but  well  damped 
oscillation  ir.  the  solution  behnet  it.  This  .ipa-lal  resolution  of  a  Jump  is 
generally  perfectly  acceptable  in  practice  (in  reality  due  to  effects  missed 
out  of  the  basic  equations  the  Jump  is  not  a  perfect  mathematical  discontinuity 
but  is  spread  over  a  short  distance),  but  if  further  resolution  is  required  an 
artificial  viscous  force  can  be  used  with  this  method  as  well.  One  aim  of  the 
present  investigation  is  to  see  what  profiles  are  obtained  for  bores  and 
hydraulic  Jumps. 

Thus,  in  summary,  the  Lax-Wendroff  method  applied  to  the  Eulerian  form 
of  the  equations  (the  method  is  also  applicable  to  the  Lagr&nglan  formulation) 
is  choson  as  the  most  promising  method  for  solving  the  shallow  water  equations; 
with  the  reservation  that  at  boundary  points  it  is  sometimes  an  advantage  co 
use  a  characteristic  method. 

Five  specific  problems  acre  chosen  for  illustrating  and  checking  the 
method,  and  for  exhibiting  some  of  the  eccentricities  of  this  type  of  flow: 

(1)  An  introductory  problem  of  a  steady  flow  down  a  channel  of 
decreasing  gradient  containing  a  hydraulic  Jump. 

(2)  A  channel  having  initially  a  steady  uniform  super-critical  flu  '. 

At  time  t  *  0  an  oscillatory  boundary  condition  is  applied  upstream.  This 
case  provides  a  check  against  a  theoretical  result. 

(3)  A  gradual  but  permanent  change  of  level  at  the  upstream  end  of  a 
channel.  The  initial  flow  can  be  either  sub-  or  super-critical,  if  the  given 
rise  of  upstream  level  is  sufficiently  abrupt  a  bore  forms  in  the  flow. 

(4)  A  similar  problem  for  a  finite  duration  surge  passing  the  upstream 
end  of  the  channel.  One  example  given  is  of  a  flood  ware  travelling  down  & 
river. 
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A  simple  hydrology  problem  with  heavy  rain  falling  on 


a  cuxiauc  vi 


variable  slope. 


The  channel  or  river  in  these  examples  is  of  uniform  rectangular  cross- 
section  and  of  width  much  greater  than  the  depth.  But  the  method  extends  in 
a  straightforward  way  to  natural  channels  in  which  the  cross-section  is 
Irregular  both  as  regard?  longitudinal  and  lateral  variation,  if  it  is  assumed 
that  the  flow  is  still  one-dimensional.  Also  the  method  is  applicable  to  one¬ 
dimensional  tidal  calculations  in  estuaries  and  channels .  The  major  effort  in 
applying  the  method  to  a  real  situation  lies  in  extracting  the  geometrical  data 
of  the  channel  from  maps  or  a  special  survey,  and  putting  it  in  a  form  suitable 
Ibr  the  computer. 

2  THEORY 

2.1  Basic  equations 

The  continuity  equation  for  winter  flowing  in  a  wide  channel  of  uniform 
rectangular  cross-section  is 


dh 

St 


+  u 


q 


0) 


and  the  momentum  equation  is 

«  0  .  (2) 

The  motion  is  assumed  to  be  one-dimensional  with  water  depth  h(x,t)  and 
velocity  u(x,t),  vhei*e  x  is  the  distance  coordinate  measured  in  the  down¬ 
stream  direction  and  t  denotes  the  time,.  Disturbances  to  the  steady  basic 
flow  are  assumed  to  produce  wuves  of  length  and  radius  of  curvature  of  the 
water  surface  much  greater  than  the  water  depth,  so  that  the  vertical 
acceleration  of  the  water  is  small  and  the  pressure  can  be  taken  at  the  hydro¬ 
static  value.  The  downward  slope  of  the  cnannel  bed  is  denoted  by  3{x)  so 
the  corresponding  slope  of  the  'water  surface  is  -  S,  giving  the  pressure 
gradient  term  included  in  (2).  The  frictional  resistance  ia  represented  by  the 
Chdzy  approximation  and,  since  the  breadth  of  the  channel  (b)  is  assumed  much 
greater  t.hp.i  the  depth,  the  hydraulic  mean  depth  (or  hydraulic  radius)  is  given 
by 


cross-sectional  area  of  water 
wetted  perimeter 


bh 
+  2h 


h  , 
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accounting  for  this  factor  in  the  denominator  of  zLx.  friction  term  of  (2).  The 
C'nezy  constant,  C,  may  vary  with  x,  to  allow  iur  changes  of  roughness.  The 
quantity  q,  occurring  in  both  (l )  and  (2),  is  the  lateral  inflow  from  rain, 
run-off  and  tributary  flow,  less  the  outflow  from  seepage,  etc.  It  is  assumed 
in  (2)  that  all  inflow  enters  the  channel  with  zero  velocity  component  in  the 
direction  of  the  main  stream,  and  that,  for  the  outflow,  this  same  velocity 
component  is  reduced  to  zero  on  leaving  the  main  stream.  The  last  term  of  (2) 
is  usually  relatively  very  small,  so  these  assumptions  are  not  crucial;  the 
dominant  effect  of  inflow  and  outflow  is  to  add  the  term  on  the  right-hand  side 
of  (l ) .  The  value  of  q  may  vary  with  both  x  and  t,  and  is  expressed  in 
units  of 

(volume/unit  time)/((unlt  length  of  channel)  x  (unit  width))  ,  e.g.  ft/s  . 

p 

Equations  equivalent  to  (l)  and  (2)  are  derived  in  detail  by  Stoker  . 

(Stoker  usee  the  Manning  approximation  for  the  frictional  resistance  which 
4/3 

introduces  h  '  *  in  place  of  h  in  the  resistance  term,  but  only  trivial 
changes  are  necessary  below  to  incorporate  this  alternative . ) 

Xf  the  bed  of  the  channel  has  constant  slope  and  the  net  inflow  term  is 
neglected,  it  can  be  seen  that  a  r--1  eady  uniform  flow  satisfies  the  simple 
relation 


u  =  C(hS)^  .  (5) 

The  simplicity  of  (3)  results  partly  from  writing  the  dimensionless  constant  in 

2 

the  resistance  term  as  g/c  . 

2 .  Equations  in  characteristic  form 

The  numerical  solution  of  ('1 )  and  (2)  is  to  be  obtained  by  a  finite 
difference  method  rather  than  a  characteristic  method,  but  it  is  useful  to 
express  these  equations  in  characteristic  form.  In  one  example  the  characteristic 
properties  are  used  directly  to  determine  an  unspecified  boundary  value,  and  a 
knowledge  of  the  slopes  of  vhe  characteristics  in  the  x,t  plane  fixes  the 
number  of  conditions  to  be  given  at  boundaries  and  also  determines  the  maximum 
value  of  the  finite  difference  time  step  for  a  given  space  step  if  the 
calculations  are  to  be  numerically  stable. 

In  the  x,t  plane  the  characteristics  of  (1)  and  (2)  have  slopes 
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^  *  u  +  (gh)*  and  ||  »  u  -  (gh)*  ,  (4) 

and  the  corresponding  characteristic  relations  are  found  in  the  usual  way 
(e.g.  Courant-5)  to  be  respectively 

£<«±W;  -  8 {s - ct;  1  ' l}  * 

»  Z  ±  ,  say  .  (5) 


Small  amplitude  disturbances  travel  both  downstream  and  upstream  at  a  wave  speed 
(gh)^  relative  to  the  water  velocity  u.  So,  for  a  Frond*  number  greater  than 

one,  i.e.  F  *»  ■  -  -r  >1,  such  disturbances  can  only  travel  downstream. 

(dO* 

3  METHOD  OF  BOLOTION 
3.1  General  method 

■n  <i  — —T» 

The  numerical  solution  of  (l )  and  (2)  is  obtained  by  a  simple  extension  of 
the  two-step  Lax-Wendroff  method.  These  equations  are  written  as 


dh 

5F 


+ 


rs 


q 


(6) 


and 

if +  k  +  ^  “  gS  ~  ^ 

or  for  brevity  as 


•nd 


where 


dh 

5t 


q 


du  dE 
5t  dx 


V  , 


0.  «  hu 


E 


+  gh 


(7) 


(8) 


(9) 


(10) 


and  y  is  the  right-hand  side  of  (7). 


A  rectangular  finite  difference  net  is  taken  in  the  x,t  plane,  with 
pivotal  points  =  i  Ax  in  the  x-direction  together  vith  a  suitable  time 
step  At,  as  depicted  in  Fig.1 .  It  is  assumed  that  the  solution  is  known  at 
all  at  time  t,  in  particular  the  initial  conditions  fix  the  solution 
at  all  x^  at  t  =  0.  The  numerical  solution  is  required  at  time  t  +  At. 

The  first  step  is  to  obtain  provisional  values  at  the  centres  of  the 
rectangular  meshes  (i.e.  points  aid-way  between  the  x^  at  time  t  +  \  At) 
from  the  explicit  formulae 

W  +  i  -  i  +  Vi  .  Vi *  «i<*>  ,  , .  - 

- JTSE - + - S -  "  *  lVi(t)  + 

....  (li) 


and 


u1+i(t  +  ?  At)  -  £  {ut(t)  +  ulx1 (t)}  Ei+1 (t)  -  Ei(t) 

- - - - ; - “+  - - 25 - 


\  (vi+1(t)  +  vi(t)} 
....  (12) 


The  second  step  uses  the^j  staggered  values  to  obtain  the  required  solution  at 
the  pivotal  points  at  time  t  +  At  from  the  explicit  formulae 

h^t  +  At)  -  h^t)  Qi+i(t  +  £  At)  -  ^(t  +  \  At) 

»  £  +  \  A*)  +  q^t  +  \  At)}  (15) 


anu 


ut(t  +  At)  -  ui(t)  B1+i(t  +  \  At)  -  E1j..(t  +  \  At) 
—  -T  Si 


4  lvi4i(t  +  4  tst)  +  Yj_j(t  +  4  At)i 


(14) 


The  finite  difference  solution  can  then  be  found  at  time  t  +  2  At,  and  so  on 
through  as  many  time  steps  as  required. 


This  method  has  the  usual  condition  for  numerical  stability: 
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It  Is  assumed  at  this  stage  that  u  St  0.  Condition  (15)  is  inferred  from  the 

•I 

corresponding  condition  for  the  equations  of  gas  dynamics  ,  and  no  instability 
has  occurred  in  the  hydraulics  examples  when  the  ratio  AxjAt  is  in 
conformity  with  (15)  at  all  points. 

3.2  Boundary  points 

In  the  examples  of  section  5,  the  river  or  channel  is  taken  to  be  in  a 
state  of  steady  uniform  flow  at  t  =  0,  with  h  =  h#  and  u  -  u„  say,  and  with 
these  constant  values  related  through  (3).  For  t  >  0,  a  disturbance  is  applied 
at  x  c  o  and  the  solution  is  soug.  downstream  of  this  point  at  subsequent 
times . 

If  F  <  1  at  x  *  0,  just  one  boundary  condition  is  required  at  this 
point.  The  explanation  being  that  only  one  of  the  characteristics  leaves  the 
region  x  >  0  of  the  x,t  plane  when  drawn  backwards  in  time,  this  charac¬ 
teristic  is  the  one  with 

|f  «  u  +  (gh)^  >  0  .  (16) 

The  other  characteristic  with 

if  »  u  -  (gh)*  <  0  (17) 

enters  the  region  x  >  0  and  intersects  time  t  at  the  positive  value  x  =  x*, 
as  shown  in  Fig. 2.  This  second  characteristic  gives  a  relation  connecting 
values  on  the  boundary  x  ■  0  with  internal  values  that  ha.  e  already  been 
found,  thus  only  one  of  the  two  boundary  values  remains  free  to  be  specified 
as  data.  Suppose  the  depth  variation  at  x  =  0  is  specified,  then  the  water 
velocity  uQ(i  +  At)  is  calculated  by  an  iterative  method  based  on  (i7)  and  the 
corresponding  characteristic  relation.  Referring  to  Fig. 2,  the  procedure  is: 

(a)  Assume  as  starting  values 

uQ(t+At)  ®  u0(t)  ,  u»(t)  =  uQ(t)  and  h*(t)  -  hQ(t) 

The  values  of  hQ(t  +  At)  and  hQ(t)  are  known  from  the  given  boundary 
condition  and  uQ(t)  hoc  already  been  found. 

(b)  Find  x3  from  ('17)  expressed  in  the  finite  difference  form 


10 


x*  =  -  f  {uQ  -  (ghQ)*  +  u*  -  (gh*)^}  At  ,  (18) 

where  u  stand s  for  u  (t  +  At)  and  similarly  for  h  . 

o  o  o 

(c)  Obtain  new  values  of  u*  and  h*  by  quadratic  interpolation 
between  the  known  values  of  u  and  h  at  (0,  t),  (x^, t)  and  (xg, t) . 

(d)  Obtain  a  new  estimate  of  UQ(t  +  At)  from  the  following  finite 
difference  approximation  of  (5)  (lower  signs  for  this  characteristic): 

uQ  =  2(ghQ)^  +  u*  -  2(gh*)^  +  \  {z_(0,t  +  At)  +  Z_  (x‘,t)}  At  .  (19) 


(e)  Repeat  from  step  (b)  on  until  the  value  of  uQ  has  converged. 


This  simple  iterative  procedure  is  not  always  convergent,  since  if  steps  (b)  to 
(d)  are  represented  by  u^  =  f  (u^)  where  u^  is  the  nth  iterate  of  uq  and 
u“+1  the  n+1-th,  then  jf*j  >1  when  At  is  sufficiently  large,  mainly  through 
the  effect  of  the  resistance  term.  In  gas  dynamics  this  term  is  absent  and 

4 

convergence  is  easier  to  obtain.  The  powerful  Wegstein  method  (Buckingham  , 
for  example)  is  used  here  to  obtain  convergence.  This  method  forms  two 

sequences  and  u”,  such  that 


u 


n 

u  - 
o 


(u"-un-1)(un-un-1) 

K  o  o  o  o' 

i  n  n-1  ”n  »1  u-2\ 

(u“-<  -u  +u  ) 


and  the  iterative  process  is  now 


f 


This  method  has  been  found  to  produce  a  convergent  value  of  uQ  (in  practice 
a  value  is  accepted  when  differing  by  less  than  0.0001  ft/s  from  the  previous 

*  —jr  *T 

iterate)  after  at  most  7  iterations.  The  starting  values  u*,  u^  at 
t  *  At  are  taken  equal  to  uQ(t). 

The  general  finite  difference  method  can  then  be  applied  to  find  the 
internal  values  h^,  itj,  etc. 

If  F  >  1  x  °  0,  both  characteristics  have  |^  >  0  and  intersect 
time  t  outside  the  region  at  negative  values  of  x.  Thus  there  are  now  no 
relations  connecting  boundary  values  with  known  internal  values,  and  both 
variables  have  to  be  given  as  data  on  the  boundary. 
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It  is  also  necessary  to  specify  the  downstream  position  at  which  the 
calculations  are  to  stop.  For  a  disturbance  starting  at  x  «  0  at  time  t  =  0, 
the  leading  disturbance  or  wavefront  travels  downstream  at  the  known  speed 

f£  -  u#  +  (gh#)*  .  (20) 

Hence,  at  each  pivotal  value  of  t,  the  calculation  can  end  at  the  highest 
value  of  i  such  that 

i  Ax  <  {u,  +  (gh*)^}  t  ,  (21 ) 

since  the  water  is  undisturbed  from  the  initial  state  at  pivotal  points  that  lie 
further  downstream. 

There  is  an  important  exception  to  this;  if  the  water  depth  at  x  *  0 
increases  sufficiently  rapidly  or  the  initial  Proud#  number  is  sufficiently 
high,  a  bore  will  eventually  fora  in  the  flow.  If  the  bore  is  at  the  head  of 
the  disturbance,  it  travela  faster  than  the  speed  given  by  (20)  and  the 
calculation  must  in  this  case  be  continued  to  a  greater  value  of  i  than 
required  by  (21 ) .  In  practice,  when  a  bore  is  likely,  a  generous  number  of 
extra  points  is  included  ahead  of  the  nornal  wavefront.  The  number  of  extra 
points  is  known  to  be  sufficient  when  the  computed  solution  at  the  last  few 
points  reproduces  the  initial  undisturbed  state. 

Further  discussion  of  initial  and  boundary  conditions  is  included  in  the 
following  examples .  In  particular,  simpler  forms  of  upstream  sad  downstream 
conditions  occur  in  sections  4  and  6. 

3.3  Solution  near  jumps 

Across  a  bore  or  hydraulic  jump  the  water  depth  and  velocity  are 
discontinuous,  and  it  can  be  shown  that  water  crossing  the  discontinuity  suffers 

C 

a  sudden  loss  of  mechanical,  energy  (lanrtr).  The  lost  energy  is  turned  into  heat 
by  turbulence  at  the  juap  (or  part  of  it  can  be  radiated  away  in  the  form  of 
surface  waves)  and  the  mechanism  for  this  abrupt  energy  loss  is  not  contained  in 
the  basic  equations. 

As  in  gas  dynamics,  introducing  an  artificial  viscous  force  term  of  the 


form 


4{(ftAx)2(S)  }' 


with  a 


constant 


f 
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gives  a  loss  of  mechanical  energy  which  is  large  where  u  is  changing  rapidly. 
The  addition  of  this  fictitious  term  enables  a  continuous  solution  to  be 
obtained  across  a  jump .  In  the  numerical  solution  the  jump  is  spread  smoothly 
over  a  short  length  of  x,  though  with  some  oscillations  behind  it,  in  place 
of  an  exact  mathematical  discontinuity.  For  a  given  value  of  Ax  the  length 
of  the  transition  and  the  amplitude  of  the  oscillations  depend  on  the  value  of 
a  that  is  chosen,  and  the  optimum  choice  is  a  matter  for  numerical  experiment. 
In  practice  this  extra  term  is  only  important  in  the  neighbourhood  of  the  jump. 

In  the  Lax-Wendroff  method  a  similar  extra  force  tern  is  provided  auto- 
matically  by  the  truncation  error,  this  being  proportional  to  d  u/dx  .  The 
distinction  between  the  differential  equations  and  their  finite  difference 
representations  on  the  Lax-Wendroff  scheme  is  very  useful,  as  it  allows  flows 
containing  jumps  to  be  covered  by  the  basic  numerical  method  without  any  special 
treatment.  This  type  of  difference  scheme  is  called  •dissipative*. 

4  STEADY  FI/3W  CONTAINING  A  HYDRAULIC  JUMP 

The  steady  forms  of  (1 )  and  (2)  with  q  -  0  are 

uh  -  Q  (constant.)  (22) 

and 

0  !  (25) 

eliminating  u  gives 


Consider  the  integration  of  (24)  in  the  downstream  direction  on  a  long 

slope  with  a  steadily  decreasing  gradient,  starting  with  a  super-critical  flow 

at  x  =  0.  As  the  gradient  decreases,  the  water  velocity  decreases  and  the 

2  3 

depth  increases,  until  a  point  is  reached  at  which  Q  =  gh  ,  that  is  F  -  1, 
and  ah/dx  tends  to  positive  infinity.  The  physically  acceptable  solution  in 
these  circumstances  consists  of  a  transition  via  a  hydraulic  jump  between  two 
branches  of  the  solution. 


WW1JL". 
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The  gradient  is  1/50  upstream  of  x  =  10  ft,  and  after  this  point  the  gradient 

decreases  smoothly  to  a  final  value  of  1  /800  downstream  of  x  =  3C  ft.  The 

corresponding  variation  of  the  depth  is  to  be  calculated.  It  is  assumed  that 
the  flow  is  unifora  in  x  SO  and  in  x  fc  60  ft  with  (3)  holding,  Q.  —  4  ft^/s, 

and  C  =  So  ft^/s.  These  conditions  fix  the  flew  at  x  =  0  as: 

h  «=  0.5  ft  ,  u  =  8  ft/s  and  F  -  2  j 

and  at  x  «  60  ft  as: 

h  *=  1.26  ft,  u  =  3.17  ft/s  and  F  =  0.5 

The  solution  has  been  obtained  by  the  artificial  viscous  force  method,  but 
for  brevity  the  details  are  not  given  here.  We  concentrate  on  the  Lax-Wendroff 
method .  The  unsteady  formulation  is  retained  and  t\  rough  guess  at  the  solution 
provides  the  initial  conditions  required  at  t  =  0.  The  time  t  now  corresponds 
to  an  iteration  parameter  and  At  is  chosen  at  each  irfcage  to  keep  well  within 
the  stability  criterion  (15)  at  all  points.  Two  boundary  conditions  are 
required  at  x  »  0  as  F  >  1  there,  these  are  h  =  0.5  ft  and  u  =  8  ft/s. 

At  x  c  60  ft,  F  <  1  so  only  one  condition  may  be  specif i  ed,  this  is  taken  as 
h  =  1 .26  ft  and  the  value  of  u  at  this  point  is  obtained  by  a  characteristic 
calculation,  similar  to  that  ot"  section  3.2.  As  the  solution  converges  to  the 
required  steady  state,  this  calculation  yJ.vee  u  -*  3. 17  ft/s  at  x  =  60  ft  as 
required. 

The  result  with  Ax  «=  1  f  t  is  shown  in  Fig. 4  after  240  iterations,  that 
is  after  240  steps  of  varying  At.  The  solution  is  hardly  changing  with  t 
at  this  stage  and  shows  the  sharp  profile  for  the  hydraulic  jump  that  is  pro¬ 
duced  by  the  method,  with  the  typical  initial  overshoot.  Also  plotted  on  Fig. 4 
ii>  the  assumed  value  of  h  at  t  =  0,  this  Is  obtained  from  h  =  Q/u  with 
u  assumed  to  vary  linearly  between  x  =  1  f t  and  x  =  59  ft.  The  values  of  u 
at  thepe  points  are  taken  equal  to  the  respective  boundary  values,  given  above. 

5  THE  PROPAGATION  OF  SISTU5BAHCFS  DCHK  CHANNELS 

5.1  Oscillatory  waves  superimposed  on  a  super-critical  flow 

We  consider  a  wide  channel  of  uniform  rectangular  cross-section  with 
constant  gradient,  3,  and  with  no  lateral  inflow.  Initially  the  flow  is  taken 
to  be  steady  and  uniform  with  depth  h#  and  velocity  u*,  such  that 
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u*  «  C(h#  S)^,  in  accord  with  (3).  The  flow  is  assumed  to  be  super-critical 
with  u#  >  (gh#)*>  so  that  S  >  g/c~.  Two  boundary  conditions  are  required  at 
the  upstream  position  x  «*  0,  these  are  taken  as 

h  *=  h#  +  H  sin  wb  ,  u  *=  u#  +  U  sin  wb  (25) 


for  t  >  0.  No  special  downstream  conditions  are  required,  but  beyond  the  wave- 
front  of  the  disturbance:  h  »  h*,  u  =  u-.  It  is  assumed  that  H/h*  and  u/u# 
are  sufficiently  small  or  F#  =>  u*/(gh*)*  is  sufficiently  large  to  keep  the 
flow  at  x  *  0  super-critical  at  all  values  of  t.  The  numerical  solution  is 
required  for  the  disturbance  produced  downstream  of  x  «  0,  in  particular  for 
the  eventual  steady  oscillation  ab  any  point. 

If  H  and  U  are  relatively  small,  the  analytical  solution  of  the 
linearised  forms  of  (j )  and  (2)  provides  a  good  test  case.  Writing  for  the 
steady  oscillation 


+  H 


«  u,  +  U  e 


i(wt-kx) 


(26) 


substituting  into  (l )  and  (2),  and  retaining  only  first  powers  of  H  and  U 
gives 


(1^  -  1 )  c*  k2  +  (3ia  F*  -  2wp»)  e,  k  +  (<*£  -  21&  w)  -  0  ,  (27) 


where 


c 


# 


and  a 


gS/u* 


(28) 


Thus 


2w  F#  -  31a  F»  ±{  (4c£  -  9a2  F*)  -  4ia  «(F?  +  2)}^ 

i  =  - - - - - 

2c*(F*  -  1 ) 

For  the  special  value  F#  =  2,  equation  (29)  has  the  two  roots 


(29) 


JL  v 
3c#  *  *2 


u  -  2ia 


(30) 


The  first  root  gives  a  wave  of  constant  amplitude  travelling  downstream  with 
speed  u*  +  c*  =  3c# •  The  second  root  gives  a  damped  wave  travelling  downstream 
with  the  lower  speed  u #  -  c*  =  c#.  Hence,  if  the  constants  in  the  boundary 
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solution  for  very  snail  H/h#  should  represent  a  wave  of  constant  amplitude. 
It  is  easy  to  verify  that  the  slower  wave  Is  absent  if 

U  * 


The  data  for  this  test  of  the  numerical  method  and  computer  programme  was 
taken  as  ? 


fc* 

«  0.5  ft 

and  hence 

c# 

*  4  ft/s  $ 

*  8  ft/s 

as 

** 

-  2  ; 

c 

-  48  ftfys 

and  hence 

s 

-  1/18  ,  a  -  2/9  j 

2a 

u 

«  0.5  « 

giving 

fa) 

=  4x  rad/s  j 

H  -  0.005  ft  {«  hJlOQ)  ,  and  U  =  0.04  ft/s  from  (31)  . 

The  finite  difference  steps  were 

At  «  0.01  s  (that  is  50  steps  of  At  in  each  cycle) 

and 

x  «=  0.15  tt  (>  (u»  +  c#)  At  -  0.12  ft) 

With  these  values  the  numerical  results  showed  an  increase  of  wave 
amplitude  of  only  0.3$  at  a  point  24o  steps  of  Ax  downstream.  The  solution 
was  advanced  over  360  steps  of  At  in  the  course  of  the  calculation .  This 
result  provides  a  useful  check  on  the  programne. 

For  general  U,  H  in  the  upstream  boundary  conditions,  it  follows  from 
(29)  that  when  F#  <  2  both  waves  are  damped  in  the  direction  of  travel.  When 
F#  >  2,  the  faster  downstream  wave  actually  increases  in  amplitude  and  even¬ 
tually  the  linear  theory  ceases  to  be  valid.  The  slower  wave  is  still  damped. 
In  practice,  instability  in  the  form  of  roll  waves  is  observed  in  steady  flows 
when  F#  >2.  Boll  waves  (Stoker  )  are  a  series  of  bores,  spaced  equi distantly, 
and  with  a  frequency  equal  perhaps  to  that  of  a  small  disturbance  applied  to 


/**  -  u*\ 
VTTSTV 


H 


e*  H 

b* 


(51) 
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the  flow,  for  example  the  frequency  of  surface  waves  in  a  reservoir  at  the 
upstream  end  of  the  channel. 

Thus,  one  interesting  application  of  the  numerical  method  for  the  solution 
of  the  non-linear  equations  is  to  see  the  results  of  applying  a  disturbance  of 
very  small  amplitude  at  x  =  0  to  a  flow  with  F#  >  2 ,  To  accelerate  the 
growth  of  the  disturbance  the  high  value  F#  =  9  was  assumed.  Thi  amplitudes 
at  x  =  0  were  taken  as 


H  =  0.025  ft  (that  is  5$  of  h#)  and  U  =  0.2  ft/a  , 


with  the  above  values  of  h*,  C  and  w.  The  results  are  shown  in  Fig. 5.  The 
depth  is  plotted  as  a  function  of  x  at  a  time  t  =3.145  s  after  the  start 
of  the  di  sturbance .  The  graph  shows  the  rapid  transition  from  small  amplitude 
sinusoidal  waves  into  sharp  peaked  waves  separated  by  water  of  nearly  constant 
depth.  This  profile  has  some  similarity  to  observed  roll  waves,  but 
unfortunately  with  this  symnetrical  steepening  the  curvature  of  the  water 
surface  is  becoming  too  large  at  the  crests  for  shallow  water  theory  to  be 
valid. 

This  sub-section  has  been  concerned  with  the  steady  oscillations 
developed  well  behind  the  wavefront.  In  fact,  for  F#  >  2  the  head  of  the 
disturbance  travelling  downstream  will  be  a  bore,  and  even  for  F#  <  2  a  bore 
develops  if  the  initial  rate  of  depth  increase  (=  H«)  is  sufficiently  large 
(idghthill  and  Wfcitham^) .  The  next  two  sub-sections  include  the  numerical 
solution  near  the  wavefront,  and  the  formation  and  behaviour  of  bores . 

5.2  Transient  effects  of  a  permanent  change  of  upstream  conditions 

We  consider  as  previously  a  wide  channel  of  uniform  rectangular  cross- 
section  with  constant  gradient  of  the  bed  and  no  lateral  inflow.  The  Froude 
number  F*  =  u^/Cgh*)*  of  the  initial  steady  uniform  flow  is  arbitrary.  One 
boundary  condition  ia  repaired  at  the  upstream  position  x  =  0  when  F#  <  1, 
and  two  boundary  conditions  when  F#  >  1 .  The  conditions  at  x  =  0  are 
assumed  for  illustration  to  take  the  simple  forms: 


and 


h  =  h*  +  H*  t  for  0  <  t  <  t* 
k  =  h*  +  H*  t* 

=  h#  +  H  ,  say  ,  for  t  fe  t®  , 


>7 


representing  a  permanent  change  of  water  level  occurring  over  time  t 1 .  In 
addition  when  F#  >  l : 

u  *  u#  +  U*  t  for  0  <  t  <  t* 
and 

u  =  u*  +  U*  t* 

«=  u#  +  U  ,  say  ,  for  tit* 


For  simplicity,  it  is  assumed  at  x  =  0  that  F  remains  either  greater  than 
one  or  less  than  one  throughout  the  motion. 

The  numerical  method  of  section  3  is  directly  applicable  to  this  problem. 
If  F  <  I  at  x  «  0,  the  special  procedure  involving  a  characteristic  is  used 
to  obtain  u  at  x  =  0. 

An  analytical  result  is  available  for  this  problem:  Light hi 11  and  Whitham^ 
show  by  expanding  the  solution  near  the  wavefront  in  a  power  series  of 
x  =  t  -  x/(u*  +  c#)  with  coefficients  functions  of  t  and  the  wavefront  being 
X  =  0,  that  a  bore  forms  when  F*  <  2  if  the  initial  steepness  (>=  initial 
discontinuity  of  dh/dt)  is  such  that 

(§)„  >  T£7  t2  ‘  T*)0  +  t')  ("  K  ’  “y)  •  (52) 


In  such  cases  the  bore  forms  at  a  time  given  by 


t 


2u# 

5s$  ->*) 


(33) 


If  F#  >  2  a  bore  always  forms  in  the  disturbed  flow  for  positive  (dh/dt )Q. 
In  the  present  example 


A  formula  that  is  useful  for  checking  the  L&x-Wendroff  treatment  of  Jumps 

is  one  for  the  speed  of  travel  of  a  bore,  separating,  say,  depths  of  h#  and 

2  5 

h  (>  h#).  If  the  bore  speed  is  5*,  then  it  may  be  shown  (Stoker  ,  Lama  )  by 
expressing  the  conditions  a?  continuity  and  conservation  of  momentum  across 
the  Jump  that  the  rate  of  volume  flow  per  unit  width  across  the  Jump  (Q), 
satisfies 
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Q2  =  \  h„(h  +  h.)  (35) 

-and 

H  =  h(fc*  -  u)  *=  h#(fc*  ~  h#)  ,  ^36) 


where  u  and  u#  are  the  water  velocities  on  the  two  aides  of  the  bore.  Hence, 
the  velocity  of  the  bore  ia 


u#  + 


I 


(57) 


showing  incidentally  that 


V  >  u*  +  (gh*)^  .  (38) 

The  numerical  solution  has  been  obtained  for  the  following  examples 
(C  =  48  ft^/s  throughout ) . 

(a)  Mvert  h,  c  8  ft,  F#  =  0.25  giving  u#  «=  4  ft/s;  and  with 
h(o,t)  increasing  from  8  ft  to  1 3  ft  in  1  h.  This  is  an  idealised  exanple  of 
a  river  subject  to  a  long  duration  flood  at  an  upstream  position.  These  values 
give  K  =  0.040  and  (dh/dt)0  *=  0.0014,  so  by  (32)  no  bore  is  predicted.  The 
numerical  solution  obtained  with  Ax  *  5000  ft  and  At  =  120  s  is  shown  in  Pig. 6 
to  be  a  wave  of  nearly  constant  profile  moving  downstream  at  about  6.9  ft/s. 
Such  a  profile  is  called  a  *monoclinal  flood  wave* . 

(b)  River;  data  as  (a)  except  the  increase  of  h(0,t)  occurs  in  only 
50  s .  This  exeLflple  corresponds  roughly  to  the  rapid  opening  or  breaking  of  a 
lock  gate.  In  this  case  (dh/dt)0  =  0.10  >  K,  so  a  bore  is  predicted.  The 
results  obtained  by  the  numerical  method  are  drawn  in  Pig. 7  and  show  the  bore 
at  the  head  of  the  disturbanct  .  These  results  were  obtained  with  Ax  =  70  ft 
and  At  =  2.5  s,  and  as  a  check  repeated  with  Ax  =  35  ft  and  At  =  1 .25  s;  the 
agreement  was  satisfactory.  The  numerical  results  give  a  final  bore  speed  of 
{*  =  22.4  ft/s,  while  (37)  predicts  22.5  ft/s  using  the  smoothed  value 

h  =  9-7  ft  just  behind  the  bore.  For  comparison,  u#  +  c#  «=  20  ft/s.  Accord- 
ing  to  the  theoretical  result  (33),  the  bore  starts  to  form  at  t  ^  85  3, 
and  the  numerical  results  show  a  bore  on  the  first  profile  after  this  time  at 
t  =  102  s. 

(c)  Steep  channel;  h*  =  0.5  ft,  F#  =  1.5  giving  u *  =  6  ft/s ;  and 
with  h(0,t)  increasing  from  0.5  ft  to  1  ft  and  u(0,t)  increasing  from  6  ft/s 
to  8.5  ft/s  in  50  8.  These  values  give  K  =  0.035  and  (dh/dt)Q  =  0.010,  so  by 
(32)  no  bore  is  predicted.  The  numerical  solution  obtained  with  Ax  =  4o  ft 
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and  At  «  2.5  s  is  shown  in  Fig. 8.  The  profile  is  similar  to  that  obtained  in 
the  sub-critical  case  (a),  and  moves  downstream  at  about  10  ft/s. 

(d)  Steep  channel}  data  as  (c)  except  the  increases  of  h(0,t)  and 
u(o,  t)  occur  in  only  5  s.  In  this  case  (dh/dt)o  *=  0.10  >  K,  so  a  bore  is 
predicted.  The  numerical  solution  obtained  with  Ax  *  4  ft  and  At  =  0.25  s 
is  drawn  in  Fig. 9#  and  shows  the  bore  at  the  head  of  the  disturbance.  The 
numerical  method  gives  a  final  bore  speed  of  £*  =11.7  ft/s,  while  (37)  predicts 
11.8  ft/s  using  the  smoothed  value  h  =  0.80  ft  just  behind  the  bore.  For 
comparison,  u#  +  c*  =  10  ft/s.  According  to  (33)>  the  bore  starts  to  form  at 
t  =  10.3  8,  and  the  numerical  results  first  show  a  definite  bore  on  the  profile 
at  t  =  15  s. 

5-3  Effects  of  an  upstream  disturbance  of  finite  duration 

The  effects  of  an  upstream  boundary  condition  (or  conditions )  representing 
a  disturbance  of  finite  duration  can  be  examined  by  a  small  modification  to  the 
programme.  The  upstream  conditions  at  x  «  0  are  taken  for  illustration  as: 

h  =  h*  +  K  sin  (*t/t*)  for  0  <  t  <  t» 
and 

h  =  h*  otherwise  . 


In  addition  when  F#  >  l : 


and 


u  =  u#  +  U  sin  (rtt/t*)  for  0  <t  <  t  • 
u  =  u*  otherwise 


The  results  (32)  and  (33)  concerning  bore  formation  apply,  with 


(39) 


in  this  case. 

The  numerical  solution  has  been  obtained  for  the  following  two  farther 
examples . 

(e)  River:  data  as  example  (a),  with  H  =  5  ft  and  t*  =1  h.  These 
values  give  K  =  0.040,  as  previously,  and  (dh/dt)Q  =  0.0044,  so  by  (32)  no 
bore  is  predicted.  The  numerical  solution  obtained  with  Ax  =  5000  ft  and 
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At  «=  120  s  is  shown  in  Fig. 10.  The  peak  of  the  disturbance  travels  at  6.6  ft/s. 
This  example  resembles  a  flood  wave  travelling  down  a  ri^er  and  the  kinematic 
wave  theory  of  Lighthill  and  Whitham^  is  applicable,  with  the  wave  property 
(downstream  waves  only)  following  from  the  continuity  equation,  plus,  in  place 
of  the  full  momentum  (or  dynamic)  equation,  a  simple  steady  state  approxi- 
nation  relating  u  and  h,  such  as  (3).  The  initial  disturbance  travelling  at 
speed  u*  +  (gh#)^  (=  20  ft/s)  is  heavily  damped  and  the  main  disturbance 
travels  (on  a  linearised  form  of  the  theory)  at  speed  j|  u*  (=6  ft/s).  Example 
(a)  is  also  amenable  to  this  approximation. 

(f)  Steep  channel1!  data  as  example  (c),  with  H  =  0.5  ft,  U  =  2.5  ft/s 
and  t*  =  15  a.  These  values  give  K  -  0.035*  aa  previously,  and 
(dh/dt)Q  =  0.105  >  K,  so  a  bore  is  predicted.  The  numerical  solution  obtained 
with  Ax  =  4  ft  and  t  «  0,25  s  is  drawn  in  Fig. 11  and  shows  the  bore  at  the 
head  of  the  disturbance.  This  solution  is  in  satisfactory  agreement  with  that 
obtained  using  Ax  =  8  ft  and  At  =  0.5  s.  The  numerical  method  gives  a  final 
bore  speed  of  ft*  *=  11.2  ft/s,  while  (37)  predicts  the  same  speed  using  the 
smoothed  value  h  =  0.70  ft  just  behind  the  bore,  ^e  extrapolations  of  the 
smooth  profile  behind  the  bore  leading  to  this  value  are  shown  dotted  on  Fig.l 1 . 
In  this  example  the  wave  profile  continues  to  rise  steeply  after  the  bore. 

For  comparison,  u#  +  c#  =  10  ft/s.  According  to  (33)*  the  bore  starts  to 
form  at  t  =  9*7  a,  and  the  numerical  results  show  a  bore  to  be  forming  at 
t  =  10  s. 

The  numerical  method  is  now  well  tested  and  is  directly  applicable  to 
flows  containing  jumps.  These  test  examples  provide  strong  support  for  its 
use  in  less  idealised  problems  which  are  completely  intractable  to  theory. 

The  computer  programmes  written  for  the  problems  in  this  Report  are 
basically  similar,  and  each  example  takes  around  15  s  of  computing  time  on  an 
ICL  1907. 

6  A  SIMPLE  ICfDROLOGY  APPLICATION 

Consider  the  situation  drawn  diagrammatically  in  Fig.l 2,  in  which  very 
heavy  rain  is  falling  on  the  surface  shown,  the  rain  starting  at  time  t  =  0. 

The  boundary  conditions  at  x  =  0  and  x  =  l  are  taken  as  u  =  0,  this  one 
condition  is  sufficient  at  each  boundary  as  F  <  1  there.  This  problem  has 
some  similarity  to  that  of  rain  falling  on  d.  cambered  road  with  blocked  gutters . 
The  slope  of  the  surface  is  assumed  to  vary  as 
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The  resulting  flow  down  the  slope  is  required  as  r.  function  of  x  and  t. 

The  usual  approach  to  ‘run-off*  problems  in  hydrology  is  by  the  kinematic 
wave  approximation^  which  combines  the  continuity  equation 

sr+lr(uh)  *  5  (41) 

with  the  assumption  that  the  flow  is  quasi-steady  and  locally  uniform  with  (3) 
holding: 


u 


(42) 


in  place  of  the  full  dynamic  equation  (2).  Or  the  more  accurate  approximation 
of  (2) 


u  <=  C 


f  (s  •  *£ 


D}* 


(43) 


can  be  taken.  J^The  forme  of  (42)  and  (43)  assume  S  4  0  and  8  i  ^ 

respectively,  otherwise  u  <=  -  C(-  hS)^  ,  u  «=  -  C  -jh  -  S^jj-  , 
respectively. J  Other  approximations  of  similar  form  have  been  proposed. 
The  combination  of  (41 )  and  (42)  gives  an  equation  of  the  form 


dh  .  3  ,  dh  „  C2  h2  dS 
5t  +  2ud£  *  q  ~  “5u“  dx 


(44) 


with  u  given  by  (42).  It  is  assumed  that  C  is  independent  of  x.  On  the 
other  hand  (41 )  and  (43)  give 


2  2 

1  „  &  _  c fjr 

2  dx  2u 


C2  h2  dS 
q  "  *2 ~  dx  ' 


(45) 


with  u  given  by  (43).  Equation  (44)  requires  an  initial  condition  and  an 
upstree.  a  condition,  but  no  downstream  condition  can  be  specified  as  the  single 
family  of  characteristics  ~  ^  allows  only  downstream  waves,  travelling 

at  the  kinematic  wave  speed  ^.u.  Equation  (45)  requires  an  initial  condition 
and  both  an  upstream  and  a  downstream  boundary  condition,  as  there  is  now  a 
diffusion  effect.  The  numerical  solution  of  (44)  can  be  obtained  by  a 
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characteristic  method,  and  that  of  (45)  is  probably  best  obtained  by  an 
iterative  application  of  a  Crank-Nicolson  type  method. 

However,  in  this  work  we  can  retain  the  full  dynamic  equation  (2)  and 
solve  by  the  Lax-Wendrcff  method.  The  initial  conditions  require  care,  the 
natural  choice  h  =  u  =  0  at  t  =  0  makes  the  right-hand  side  of  (7) 
indeterminate.  One  alternative  is  to  assume  an  initial  non-zero  water  depth 
that  is  equi valent  to  &  few  seconds  of  rainfall,  with  u  given  initially  by 
the  positive  root  of  the  quadratic 

gS  -  «  0  ,  (46) 

C  h 

that  is  the  right-hand  side  of  (7)  put.  equal  to  zero.  Other  starting  assump¬ 
tions  may  be  preferable.  It  is  in  any  case  assumed  uncritically  in  this  problem 
that  the  basic  equations  (l )  and  (2)  hold  in  the  early  stages  of  the  motion 
when  the  depth  is  very  small. 

The  boundary  conditions  are  u  =  0  at  x  =  0  and  x  ®  t.  As  the  values  of 
h  at  t-hese  points  can  be  obtained  by  a  small  modification  to  the  general  finite 
difference  method,  the  characteristic  method  is  not  necessary  in  this  special 
case.  The  value  of  At  is  chosen  at  each  step  to  keep  well  within  the  stability 
requirement  (15). 

Consider  as  a  numerical  example: 

q  -  10'4  ft/s  ,  corresponding  to  about  4  inches  of  rain  falling  in 
one  hour  with  no  surface  seepage. 

I  =  100  ft,  with  a  maximum  slope  S(£A)  -  5/100  , 

C  =  48  ft^/s  and  the  initial  depth  is  taken  to  be  equivalent  to  50  s 
of  rainfall. 

The  results  obtained  with  Ax  ~  4  ft  are  shown  in  Fig. 15*  A  partial  check  can 

be  applied:  total  volume  of  rain  fallen  per  unit  length  of  surface  normal  to  the 

2 

curve  of  Fig. 12  at  t  -  116  s  equals  (116  +  30)  X  10  X  100  =  1  .46  ft  j  and 
in  fact  the  area  under  the  corresponding  curve  of  Fig. 13  agrees  with  this  figure. 
The  oscillations  t.he  solution  become  relatively  less  as  t  increases,  and 
may  be  partly  due  to  the  artificial  initial  condition  of  uniform  depth  which 
omits  a  hydraulic  jump,  chough  implying  a  deceleration  from  super-critical  to 
sub-critical  flow  at  about  x  =  So  ft.  The  oscillations  are  reduced  by  taking 
At  much  smaller  at  the  start  of  the  calculations  than  is  required  by  (15).  It 
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found  tlxat  the  relationship  between  h  and  u  in  this  region  is  close  to  that 
of  the  kinematic  approximations  (42),  (43)  and  (46).  There  is  little 
difference  between  these  three  approximations  in  this  example.  Thus,  the 
kinematic  theory  is  adequate  in  x  <  8o  ft,  but  in  x  >  8o  ft  the  dynamic 
theory  is  necessary  for  this  rather  demanding  test  case  involving  very  rapid 
deceleration,  partly  by  a  hydraulic  jump,  as  the  flow  approaches  the  downstream 
boundary . 

7  EXTENSIONS  TO  NATURAL  CHANNELS  AND  TO  ONE-DIMENSIONAL  TIDAL  CALCULATIONS 

In  a  river  or  estuary  the  cross-section  of  the  channel  is  rarely 
rectangular  and  uniform  in  the  along-stream  direction,  x.  The  width  at  the 
water  surface,  for  example,  can  vary  with  both  x  and,  at  a  fixed  value  of  x, 
with  the  current  level  of  the  water  surface.  The  water  depth  varies  in  the 
across-streara  direction  and  often  the  variation  of  the  bed  in  the  along-stream 
direction  is  such  that  even  the  mean  water  depth  does  not  vary  smoothly  with 
x.  However,  the  height  of  the  water  surface,  n(x, t),  above  a  fixed  hori¬ 
zontal  plane  (for  example  mean  sea  level)  varies  much  more  smoothly  with  x 
and  lias  hardly  any  variation  across  the  width  of  the  stream,  and  for  these 
reasons  is  chosen  as  a  dependent  variable  in  the  general  case.  See  ?lg.l4. 

In  fact,  for  the  rectangular  channel,  rj  was  eliminated  from  (2)  through  the 
relation 


dr)  _  dh 
dx  dx 


S 


Assuming  that  the  flow  is  still  approximately  one-dimensiona  ,  (6)  and  (7) 
generalise  to 

ST +  5x  (^u)  =  q 


and 


du 

dt 


+ 


+  en) 


C  R 


-  2E 
A 


(48) 


where  A(x,t)  is  the  cross-sectional  area  of  the  water  at  location  x  and 
time  t,  aid  R(x,t)  is  similarly  the  hydraulic  mean  depth.  The  lateral 
inflow  q  is  now  in  units  of 
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These  equations  are  equivalent  to  those  derived  by  Stoker  . 

For  the  uniform  rectangular  channel,  A  was  simply  proportional  to  h. 

But  now  at  each  value  of  x  there  is  a  functional  relationship  between  A 
and  n  to  be  determined  from  large  scale  maps  or  charts,  or  perhaps  from  a 
special  survey.  This  is  depicted  diagrammatically  in  Fig.l 4,  which  shows  u 
as  a  different  function  of  A  at  each  pivotal  value  of  x.  Similarly,  R  is 
obtained  as  a  function  of  A  for  each  x. 

If  the  data  of  Fig.l 4  is  given,  the  numerical  solution  can  be  obtained  by 
the  Lax-Wendroff  method  as  before,  the  dependent  variables  now  being  A  and  u. 

At  several  stages  of  the  calculation  values  of  n  and  R  have  to  be  obtained  from 
this  data  corresponding  to  known  values  of  A  and  x. 

In  estuaries  and  in  tidal  rivers  and  channels,  the  water  velocity  may 
change  in  direction  over  the  tidal  cycle.  The  following  points  now  apply: 

2 

(a)  The  factor  u  in  the  frictional  resistance  term  must  be  written 
as  |u|  u,  to  ensure  the  resistance  force  always  opposes  the  motion.  This 
requires  only  a  simple  change  in  the  programme. 

(b)  The  stability  criterion  (for  a  rectangular  channel)  is  now 

Ax  >  (|u|  +  (i?h)^)  At  , 


in  place  of  (15). 

(c)  Usually  F  <  I  at  both  the  mouth  and  the  landward  limit  of  the 
estuary,  thus  one  boundary  condition  is  required  at  each  point.  For  example, 
the  tidal  level  might  be  specified  at  the  mouth  as 

h  =  h*  +  H  sin  wc  , 

and  the  corresponding  velocity  is  obtained  by  the  method  of  section  3.2;  while 
at  the  landward  limit  there  might  be  a  barrier  with  the  condition  u  =  0, 
and  the  corresponding  water  level  is  obtained  as  in  section  6. 

(d)  Approximate  initial  conditions  must  be  specified  at  t  =  0  and  the 
calculations  then  advanced  through  sufficient  time  steps  for  the  solution  to  have 
become  periodic  in  t.  The  better  the  guess  at  these  conditions,  the  quicker 
this  will  be  achieved,  but  very  rough  initial  conditions  may  suffice  in  practice. 
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method  covers  this  case  without  any  special  procedures  in  the  computer 
programme . 


8  CONCLUSIONS 

The  Lax-Wendroff  method  has  been  applied  in  this  Report  to  a  selection  of 
one-dimensional  unsteady  problems  of  open  channel  hydraulics.  It  has  been 
shown  that  the  method  can  be  applied  to  flows  containing  bores  and  hydraulic 
jumps  without  either  shock  fitting  or  employing  an  artificial  viscous  force, 
and  at  the  same  time  gives  an  acceptable  spacial  resolution  to  these 
discontinuities .  The  method  is  simple,  easy  to  programme  for  a  computer,  and 
even  for  flows  without  jumps  has  advantages  over  other  methods.  The  extensions 
of  section  7  cover  many  practical  situations,  and  further  extensions  are 
possible,  for  example  junctions  and  two-dimensional  unsteady  problems  could  be 
considered. 
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dimensionless  constant  in  artificial  viscous  force 
area  of  water  cross-section 
channel  breadth 

wave  speed  in  undisturbed  water:  (gh*)^ 

Chezy  constant 
|u2  +  gh 

finite  difference  approximation  to  E  at  x  =  x^ 

Froude  number:  u/(gh)^ 

1’roude  number  of  undisturbed  flow 
acceleration  due  to  gravity 
water  depth 

finite  difference  approximation  to  h  at  x  - 

undisturbed  water  depth 

value  of  h  at  x  =  x* 

amplitude  of  depth  disturbance  at  x  **  0 

rate  of  increase  of  water  depth  at  x  =  0 

square  root  of  -1  in  section  5*1 

wave  number  in  (26) 

roots  of  (27) 

defined  by  (52) 

a  fixed  downstream  boundary  is  taken  as  x  =*  l 
lateral  inflow 

finite  difference  approximation  to  q  at  x  = 

rate  of  volume  flow  per  unit  breadth:  uh 
finite  difference  approximation  to  Q  at  x  »  xi 

hydraulic  mean  depth 
downward  slope  of  channel  bed 
time 

duration  of  change  of  upstream  conditions 
water  velocity 

finite  difference  approximation  to  u  at  x  *=  xi 

undisturbed  water  velocity 
value  of  u  at  x  =  x* 
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u“  u"  the  two  sequences  used  in  the  Wegstein  method  of  section  3.2 
(n  *  0,1,2,  . . . . ) 

U  amplitude  of  velocity  disturbance  at  x  *=  0 

U*  rate  of  increase  of  water  velocity  at  x  »  0 

V  abbreviation  for  right-hand  side  of  (7) 

V^t)  finite  difference  approximation  to  V  at  x  *  x^ 

x  distance  along  the  channel  in  the  downstream  direction 

x^  finite  difference  mesh  points  in  the  x-direction  (i  *  0,1,2,  . ...) 

x*  the  characteristic  with  negative  slope  from  (O,  t  +  At)  passes 

through  (x*,t)t  see  section  3.2  and  Fig. 2 
Z.Z  abbreviations  for  right-hand  sides  of  (5) 

T  «» 

a  g^u# 

Ax,  At  finite  difference  steps 

t]  vertical  displacement  of  water  surface  from  a  given  horizontal  plane 

velocity  of  bore 

x  characteristic  variable  defined  in  section  5.2 

w  circular  frequency  of  oscillatory  disturbance 
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